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Abstract 

Polarizable interaction potentials, parametrized using ab initio electronic struc- 
ture calculations, have been used in molecular dynamics simulations to study the 
' (—^ , conduction mechanism in Y2O3- and SC2O3 -doped zirconias. The influence of 

I vacancy-vacancy and vacancy-cation interactions on the conductivity of these ma- 

C terials has been characterised. While the latter can be avoided by using dopant 

flj cations with radii which match those of Zr + (as is the case of Sc +), the former is 

r^ an intrinsic characteristic of the fluorite lattice which cannot be avoided and which 

O is shown to be responsible for the occurrence of a maximum in the conductivity 

C/3 at dopant concentrations between 8 and 13 %. The weakness of the Sc-vacancy 

O interactions in Sc203-doped zirconia suggests that this material is likely to present 

C/3 the highest conductivity achievable in zirconias. 

,Ph 1 Introduction 

^-H Yttria-doped zirconia (YSZ) is a fluorite-structured oxide ion conductor used as an 
^ electrolyte in a solid-oxide fuel cell (SOFC) inter alia. Substitution on the cation 
t sublattice of the Zr'^+ions by the larger Y^+ ions stabilizes the eight-coordinate cu- 
bic structure and, because of the lower valency, introduces oxide ion vacancies on 
the anion lattice which then allow oxide ion diffusion.'^' If diffusion occuiTed by an 
independent vacancy-hopping mechanism, one might expect that the dependence of 
the conductivity on vacancy concentration, c, would be c(l — c), which maximizes at 
( ^) c = 0.5. However, the maximum actually occurs at a much lower vacancy concentra- 
T-H tion - where roughly 4% of the oxide ion lattice sites are unoccupied PI This has been 
L| ascribed to the reduction in the vacancy mobility by ordering processes. These can be 
• ^H caused by the differences between the host and dopant cations, which lead to a pref- 
/\ erence for vacancies to bind close to a certain cation species,'2icJ or by the interaction 
J;S —. 
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between the vacancies themselves. '^'2] 

In the case of YSZ, adding Y^+ ions has two competing effects: on the one hand, 
adding more dopant increases the number of vacancies and therefore increases the mo- 
bihty of the remaining anions; on the other hand, since the migration barrier of an 
oxygen crossing a Y-Zr or Y-Y edge is much higher than in the case of a Zr-Zr edge™ 
(because Y^+ is a much bigger cation-" than Zr^+), adding more dopant will increase 
the fraction of cation edges associated with high migration energy barriers, thus effec- 
tively reducing the overall number of mobile oxygens. It is worth noting that this can 
be seen from a vacancy perspective: vacancies prefer to bind to be nearest neighbours 
to the smaller cations (Zr'*+ in this case) and this reduces their overall mobility. If we 
express the oxygen diffusion coefficient, D, as " 

D^Doexpi-^) (1) 

where Do is a pre-exponential factor, T the temperature in Kelvin, ks Boltzmann's con- 
stant and A//„, is the oxygen migration enthalpy, increasing the number of vacancies 
will increase Do whereas increasing the number of Y-Zr or Y-Y edges will increase 
the migration enthalpy, AH,„. The resulting diffusion coefficient will peak for a certain 
dopant concentration, which is smaller than would be expected from an independent 
vacancy-hopping mechanism. This explanation is almost universally accepted^ ^ 11^ 
and has been generalised to other doped fluorite systems.l^IIStlSl in scandia-stabilised 
zirconia, however, the two cat ions have very similar radii and cation-vacancy ordering 
effects are very wealcl^l^l^^l^ For this reason the drop in the conductivity observed in 
this material (at a slightly greater dopant concentration than in YSZ) has to be caused 
by some other effect. 

Vacancy-vacancy interactions h ave a lso been proposed as a factor hindering the 
ionic conductivity of these materials F ' ^ ^"^'^ iDirec t evidence of vacancy- vacancy order- 
ing even at the high temperatures of interest for conductivity studies, is seen in diffuse 
neutron scattering patterns obtained for YSZ. ^ In this study, intense diffuse peaks are 
observed in the neutron and x-ray diffraction patterns at positions which are forbidden 
for a simple fluorite lattice. Some of these peaks have been interpreted as caused by the 
local relaxation around small aggregates in which pairs of vacancies along the (111) 
direction pack together along the (112) direction. That this vacancy pair configura- 
tion is favourable is validated by the ab initio static energy calculations of Bogicevic 
and co-workers'^ and Pietrucci et al^ By studying the relaxation times observed in 
the quasielastic neutron scattering at different points in reciprocal space, Goff et al. 
showed that the vacancies associated with the aggregates moved more slowly than iso- 
lated vacancies and suggested that this was the cause of the decrease in the conductivity 
at high vacancy concentrations. 

In the preceding paper we examined cation-vacancy and vacancy-vacancy order- 
ing effects in the Zro.8(Y/Sc)o.20i.9 system, using computer simulations. We studied 
systems in which the proportion of Y^+ and Sc^+ is varied at constant vacancy concen- 
tration (this is 5%, close to the composition at which the conductivity maximum occurs 
in YSZ). The simulations reproduced the structural information from diffraction exper- 
iment and also the variation of the conductivity as the Y/Sc ratio was varied. We found 
evidence of both cation-vacancy and vacancy- vacancy ordering and the associated en- 



ergies were consistent with the ab-initio data.^ The cation- vacancy effects were much 
stronger in the Y-rich samples, which is consistent with the lower conductivity of those 
materials, whilst the vacancy-vacancy ordering seemed to be almost independent of the 
nature of the dopant cation, at least at these dopant levels. 

The objective of the present paper is to examine the interplay between these two 
mechanisms and to discover how they combine to determine the conductivity of a par- 
ticular material at a particular doping level. We will examine the properties of (Y203).v 
- (Zr02)i-A and (Sc203),v - (Zr02)i-A (which we dub ScSZ) obtained in simulations 
at high temperature, with x varying and, with it, the vacancy concentration. The ulti- 
mate purpose is to elucidate to what extent it might be possible to further improve the 
conductivities of this class of material by further optimisation of the doping strategy. 
Although both mechanisms clearly exist, it is not clear to what extent they are cou- 
pled. Does the cation-vacancy interaction influence the vacancy-vacancy interaction 
to a significant degree, for example? Our strategy is, firstly, to compare the results of 
our simulation with measurable quantities which bear upon the properties of the va- 
cancies, notably the conductivity and diffuse scattering. No new experimental data is 
reported in the present paper, the conductivities of both systems have been measured 
previouslyl^l^ and the diffuse neutron scattering in the YSZ system was measured in 
single crystal studies by Goff et al. ^ for 0.09 < x < 0.25. Secondly, having validated 
that our simulations are reproducing these observables sufficiently closely to assert that 
the behaviour of our simulated systems is paralleling that of the real materials in this 
regard, we will examine directly the properties of the vacancies themselves in a way 
which cannot be replicated in a real experiment. In addition to the simulation of these 
materials with realistic potentials, we will also report results for model systems for 
which the potentials have been altered (for example, by equalization of cation charges) 
in order to illustrate the physical factors responsible for a particular aspect of the ob- 
served behaviour 



2 MD simulation and analysis 

The interaction potential (which we usually dub DIPPIM), and the procedure we used 
to parameterize it, has already been described in preceding work.l^^JtSSI jjj fjjjs poten- 
tial, the ionic species carry their valence charges (Zr^+, Y^+, Sc^+ and O^^), and the 
polarization effects that result from the induction of dipoles on the ions are accounted 
for. The parameters of the interaction potential for these systems were obtained from 
the application of a force-and dipole-matching procedure aimed at reproducing a large 
set of first-principles (DFT) reference data^^ on the condensed phase. Such potentials 
have been shown to provide an accurate and transferable representation of the interac- 
tions in a number of oxides. I^SSSl 

In addition to these "realistic" potentials we also construct "ideal" model systems 
(which we will refer to as i-YSZ^) which contains the same concentration of vacancies 
as (Y203)a - (Zr02)i-., with the same x value. In i-YSZ^ all the cations carry the same 
charge, such that the total cation charge balances that of the O^^ ions present in the 
simulation, and all cations have the same short-range interaction potentials as the Zr^+ 
ions in the original potential. This idealized system allows us to eliminate the effects of 
differing charges and lattice strain induced by having host and dopant cations and gives 



an idealised reference system in which the cation-vacancy ordering effects are absent. 

All the simulations were performed, unless otherwise stated, using a cubic sim- 
ulation box with 4x4x4 unit cells of the fluorite structure, i.e. 256 cations and a 
variable number of oxygen ions, depending on cation composition. For each dopant 
concentration, a certain number of Zr^+ ions were replaced with Sc^^ or Y ^+ ions and 
their positions were randomly distributed over the cation sublattice .123211 In our pre- 
vious work'^^' we concluded that local cation ordering, as observed in YsNbO?, was 
much less important in Xi^J^i^i and we presume that this remains true at other com- 
positions as well as in ScSZ (where the site mismatch between Sc^+ and Zr'*+ is even 
smaller than in the X^^ and Zr'^^ case'^' ), so that a random distribution of cations is 
appropriate. Vacancies were randomly assigned to the oxide sublattice in the initial 
configuration. The time step used was 1 fs and all the runs were performed at constant 
volume and temperature (NVT ensemble), with thermostats as described elsewhere.'^ 
The cell volume was obtained from a previous run in a NPT ensemble with zero applied 
pressure. Coulombic and dispersion interactions were summed using Ewald summa- 
tions while the short-range part of the potential was truncated to half the length of the 
simulated box (usually about 10 A). 

In order to calculate ionic conductivities, we ran long simulations (0.5 ~ 5 ns) at 
high temperatures (T > 1250 K) on YSZ and ScSZ. Each simulation was long enough 
to allow each oxygen ion to move, on average, by a distance of at least ^ 2.6 A, i.e. the 
oxygen-oxygen bond distance. Ionic conductivities were estimated from the slope of a 
plot of the mean-square displacement of the charge versus time, i.e)^ 

A = -^— limt-\\ yq,5ri{t) \^), (2) 

6kBTV t^oo \i^^' 'V 7 I /' 

where qi is the charge on ion /, 5r,(f) the displacement made by the ion in time t 
and (....) denotes an average over the simulation run. This quantity suffers from poor 
statistics, especially when the mobility of the ions is low, and so we have also calculated 
conductivities from the Nernst-Einstein expression 

A = ^^limr'£^2^|5r,(0P), (3) 

which neglects correlations between the diffusive jumps of different ions. 

To compare with the diffuse neutron diffraction results we calculate the intensity of 
(total) scattering at some point q in reciprocal space from 

/(q) = (EEa,-«>"'-'-'^); (4) 



where the sums run over all ions in the sample, r,y is the vector joining ions / and 
j and Ui is the neutron scattering length of the species to which / belongs. Because 
of the periodic boundary conditions, the accessible vectors q with a cubic simulation 
cell of side L are restricted to the set j^{m,n,p) where m, n, and p are integers, and 
this provides a limit to the resolution of the simulated pattern. For these calculations 
we used a cubic simulation box with 6x6x6 fluorite unit cells, which gives a res- 
olution of approximately 0.2 A^^. We made extensive use of the cubic symmetry to 



average over the scattering calculated from symmetry-related q-vectors. These simu- 
lations were started at 1800 K and then slowly cooled down to 800 K with a cooling 
rate of lO'^ K s^'. In this temperature range, oxygen ion diffusion is observed while, 
below 800 K, no diffusion can be observed on the available timescale. For this reason 
the simulations at low dopant content (x = 0.9) were further quenched down to room 
temperature (300 K). Although this should not influence the vacancy ordering this will 
increase the tetragonal distortions in the sample (see discussion below). 

The relaxation of the structures responsible for the diffuse scattering at a partic- 
ular point in reciprocal space q can be measured in a quasielastic neutron scattering 
experiment. In the simulation the quantity to be calculated for comparison with the 
experimental data is the intermediate scattering function /(q,f )' given by 

again, the calculated quantities may be averaged over the symmetry-related q-vectors. 

We can identify the positions of vacancies for some instantaneous ionic configura- 
tion by finding which of the coordination tetrahedra around the anion sites are empty. 
The details of this procedure are reported in previous papers .1^51221 gg^.^^gg jj^g cations 
are not diffusing we can monitor the properties of each tetrahedron from the identities 
and instantaneous positions of the four cations which sit at its vertices. Such a tetrahe- 
dron is empty if no anion is within the volume bounded by the four planes defined by 
the positions of each set of three of the cations. Because the oxide ions may perform 
large-amplitude vibrations about their average sites at the temperatures of interest for 
dynamical studies we only assign a vacancy to a tetrahedral site if the tetrahedron has 
been empty for a minimum of two frames (i.e. 100 fs). The position of the vacancy 
is defined by that of centre of the tetrahedron (which is given by the average positions 
of the four surrounding cations). Once the positions of vacancies are identified, we 
can us them to build radial distribution functions (rdfs) to study the ordering tenden- 
cies in real space. Integration of the rdf can be used to define coordination numbers. 
For example integrating the vacancy-vacancy rdf, gv-v from zero out to the position 
re of first minimum of the gv-v rdf gives the average number of vacancies which are 
nearest-neighbours another vacancy: 

«v-v=47rp/ drr gv-v{r), (6) 

Jo 

where p is the density of vacancies. 



3 Comparison with experimental data 

In this section we will compare our MD data with the available experimental data on 
(Y203)v-(Zr02)i-x and (Sc203)x-(Zr02)i-j:. Although our methodology has already 
been thoroughly tested in the preceding papep^^ as well as in our previous work^^ 25] 
for jc = 0.11, jc = 0.33, we will first ensure that our simulations at different dopant 
concentrations reproduce the experimental data which is influenced by the vacancy 
ordering effects, before drawing conclusions on the conduction mechanism of YSZ 
and ScSZ. A comparison will be therefore made with the conductivity, quasielastic and 



diffuse scattering data.l^SlESlSl 

3.1 Conductivity 

In figure [T] we show the calculated Nerrnst-Einstein conductivities for both YSZ and 
ScSZ as a function of dopant concentration at 1670 K. For YSZ, we also report the 
experimental values from ref.!^^!^ The MD data on YSZ are seen to be in reasonable 
quantitative agreement with experiment. The dependence of the conductivity on the 
dopant concentration also looks promising. Both curves show a peak in the conduc- 
tivity at jc = 8-9 % and 12-13 % for YSZ and ScSZ respectively in good accord with 
the experimental data. However, the rise in conductivity shown by our simulations is 
smaller than that of the experimental data. This may be caused by two effects. Firstly, 
at higher dopant concentrations, the Nernst-Einstein approximation of independent ion 
jumps is not accurate. To this end, we calculated the ionic conductivities from equa- 
tion l2] at the dopant concentrations shown in figure [T] and found that these are in close 
agreement with the Nernst-Einstein values for ;c < 0.10 but systematically lower to an 
increasing degree as x increases, forx >0.10. Secondly, the small systems studied in 
this paper, which do not allow the formation of grain boundaries and tetragonal do- 
mains, might affect the behaviour of the low dopant concentration part of the curve 
(see also discussion below). 

The conductivities shown above were calculated from simulations in which the cu- 
bic fluorite symmetry was enforced. This was done by equilibrating the simulation 
using an isotropic barostat,^^ before conducting the production runs at constant vol- 
ume. As it is well known, ' ^ YSZ and ScSZ are not cubic at low temperatures and 
dopant concentrations. For this reason, when the conductivity versus dopant concen- 
tration curves are calculated at lower temperatures, a slightly different behaviour is 
observed (see bottom part of figure [T]i. The maximum has, in fact, shifted to lower 
dopant concentrations (approximately 6 % and 8 % for YSZ and ScSZ respectively) 
and the curve's shape looks very similar to that of Gd203-doped Zr02 (see figure 4 
in reference, a material which, due to Gd's big ionic radius, is cubic for almost all 
compositions). A shift o f the m aximum's position as the temperature is lowered is also 
observed experimentally,'22i34| though ^^^ ^^ strong as that indicated by the comparison 
of the top and bottom part of figure [1] 

We think that the changes observed in figurefllare mainly caused by the fact that the 
real material is not locally cubic at these temperatures and low dopant concentrations^^ 
, so that the simulations, in which the cubic symmetry is enforced within a cell which 
is too small to allow the formation of distorted domains, tend to overestimate the con- 
ductivity of these materials. To confirm this, we relaxed the cubic fluorite symmetry in 
a simulation of 4% ScSZ and allowed the simulation cell to change its shape, by using 
an anisotropic barostat.l^ The simulation quickly adopted a non-cubic symmetry. The 
conductivity obtained for this system is shown in the bottom part of figure[T]as a green 
diamond. It can be appreciated that this is indeed much lower (more than ten times!) 
than the conductivity obtained from a cubic simulation. In the reminder of this paper, 
however, we will focus on highly-doped stabilised zirconias only, i.e. those materials 
with more than 8 % of dopant, which are fully stabilised at every temperature. These 
systems are therefore not affected by the above-mentioned problem and therefore sim- 
ulations, in which the cubic symmetry is enforced, can be used. 
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Figure 1 : Top: MD conductivities vs doping (x) for YSZ (red curve) and ScSZ (black curve) 
at 1670 K. We also report the experimental data frorr i i^ 'for YSZ at the same temperature. 
Bottom: MD conductivities vs doping (x)for YSZ (red curve) and ScSZ (black curve) at 1250 K. 



3.2 Diffuse scattering in YSZ 

As we remarked in the Introduction, the most direct evidence for vacancy-vacancy or- 
dering effects in the temperature range of interest for conductivity measurements comes 
from the single-crystal diffuse neutron scattering data obtained by Goff ef al. ^ In figure 
[2]we have compared the experimental patterns obtained for (¥203)^ - (Zr02)i-.v> with 
X = 0.09 and x = 0.24 respectively, with patterns calculated in the simulations for very 
similar compositions. 
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Figure 2; Experimental (left) and computed (right) neutron diffraction patterns for 9 % (top) 
and 24 % (bottom) YSZ. 



The X = 0.09 (T = 300 K) composition corresponds to the maximum of the conduc- 
tivity plot for YSZ. The figure shows a (110) plane in reciprocal space. Sharp intense 
peaks are seen in the experimental data at the fee Bragg positions (002), (004), (111), 
(113), (220), (222), etc.; in the simulation data the Bragg peaks are centred on the 
same q values but are broader because of the lower resolution. The remainder of the 
intensity which appears in the pattern is diffuse scattering and is indicative of local 
structural deformations. At x = 0.09 the most intense diffuse features in the experi- 



mental data appear at (114) and (112), these are forbidden reflections for the fluorite 
structure. According to the analysis of Goff et a. I these peaks are a consequence of 
local tetragonal distortions of the fluorite structure which occur at relatively low con- 
centrations of Y2O3. Recall that in pure Zr02, a tetragonal phase is more stable than 
the cubic one and that Y2O3 is added to stabilize the latter. At low Y2O3 concentra- 
tions local tetragonal distortions with random orientation occur even within a single 
crystal. Comparison with the MD data shows that the same diffuse features are present 
in the simulation data, together with other features which are comparatively weak in 
the experimental data. This indicates that the tendency to local tetragonal distortion is 
reproduced in the simulated system: the differences in the intensity distribution and the 
widths of the diffuse features compared to experiment may be attributed to the small 
size of our simulation cell which sets a limit to the range of correlation of the tetrago- 
nal deformation which can be accommodated and thus leads to broader, weaker diffuse 
features. Both sets of data also show a broad, weak feature at q = (1.7, 1.7, 1) which is 
associated with a pattern of lattice deformation about isolated vacancies. ^ 

As the Y2O3 concentration is increased to x = 0.24 (T = 800 K), the (114) and 
(112) features disappear and the diffuse scattering becomes dominated by new fea- 
tures (highlighted by the arrows in the bottom part of figurel2]i which may be described 
as belonging to a superlattice at G± (0.4, 0.4, ±0.8), where G is an fee Bragg peak po- 
sition.'' This pattern has been interpreted as caused by small aggregates of divacancy 
clusters (vacancies paired along the (111) direction) packed along the (112) direction, 
as found in the Zr3Y40i2 compound. These aggregates are typically 15 A in diam- 
eter.!^ This structure is consistent with the most stable arrangement of vacancy pairs 
found by Bogicevic ^ in his analysis of vacancy-ordering tendencies in YSZ and with 
our analysis in the preceding paper. Comparison with the simulation data shows that 
the signature of this particular vacancy-vacancy ordering has been reproduced in the 
simulation. Interestingly, since our simulations at x = 0.24 were run at T = 800 K, this 
means that these ordering effects are still present at high temperature and that therefore 
they will affect the conducting properties of these materials at the temperatures of in- 
terest for technological applications. 

Although we have only discussed the diffuse scattering results for YSZ itself, sim- 
ilar diffuse scattering patterns were obtained in calculations on ScSZ, for which there 
is no single crystal data. 

3.3 Quasielastic scattering 

As a final comparison with experiment, we now attempt to reproduce the quasielas- 
tic data from ref. Goff et al. showed that the relaxation times of the quasielastic 
scattering at different points in reciprocal space varied significantly. Close to the dif- 
fuse scattering feature associated with the isolated vacancies (at q= (1.67, 1.67, 1) 
the intermediate scattering function at 1800 K relaxed quite rapidly, on a timescale 
of 1.3 ps, which gives a timescale for the independent vacancy hopping. However, at 
q= (1.33,1.33,1.67) and q= (1.33,1.33,1.67), where the diffuse scattering is associated 
with clusters of vacancy pairs, the relaxation was much slower, leading to an unre- 
solved elastic peak in the associated spectra. Goff et al. deduced from this that those 
vacancies involved in the clusters moved much more slowly than the free ones, giving 
compelUng evidence for the effect of vacancy -vacancy interactions on diffusion. 
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Figure 3; Time scan at constant q through the quasielastic diffuse scattering from 13 % YSZ at 
1800 K. 



The intermediate scattering function at these three q values was calculated as ex- 
plained in the previous section, from a simulation at 1800 K on 13 % YSZ; these were 
the conditions of the experimental study. The results are shown in figure |3] It can be 
seen that only at q = (1.67, 1.67, 1) does the intermediate scattering function relax to 
zero on the timescale of the simulation. At q = (1.67, 1.67, 1) the relaxation time is 
1.7 ps^ to be compared with the 1.3 ps seen experimentally. At the other q values 
there would be a significant unresolved "elastic" peak in the spectrum indicating that 
the defect aggregates responsible for the scattering at these q values are immobile on 
the timescale on which the correlation function has been calculated. The simulations 
are thus in good quantitative accord with the experimental results. 

4 Analysis of the vacancy ordering 

Although the diffuse scattering comparison discussed above demonstrates that the or- 
dering effects present in the real material are reproduced in the simulations, the degree 
of ordering is more conveniently illustrated in real-space, by comparing cation- vacancy 
and vacancy-vacancy radial distribution function, obtained for different dopant concen- 
trations. 

In the preceding paper, we have already illustrated how the cation-vacancy radial 
distribution functions provide information about the tendency of vacancies to associate 
with particular cation types. Our results are consistent with earlier suggestions'^n^'^^l 
that the cation- vacancy interactions are driven by the lattice strain associated with the 
difference in size between the host and dopant cations. The fact that vacancies and 



llxhis was obtained by fitting an exponential decay to the curve in figure 3 
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trivalent dopants have the same charge, and would therefore be expected to repel, seems 
to play a much smaller role. It was shown that in YSZ, where there is a substantial dis- 
parity in size between the two cations, there was a high tendency for the vacancies to 
occur in the first coordination shell of Zr, relative to Y. This can be interpreted as allow- 
ing the Zr to reduce its coordination number from the 8 of the normal fluorite structure. 
In ScSZ, where the cation sizes are much closer, this tendency was much less marked. 
The stronger cation-vacancy ordering effect in YSZ provides a natural way of explain- 
ing the higher conductivity of ScSZ at the low dopant level which was the focus of 
the previous study. However, in view of the weakness of the effect in ScSZ, unless 
the cation-vacancy becomes much stronger at higher dopant levels, it seems likely that 
other factors contribute to the drop in conductivity and the appearance of a maximum 
at relatively low dopant levels. 

In figure |4] we show the cation-vacancy radial distribution functions for YSZ and 
ScSZ for different dopant concentrations at 1250 K. We choose to show the concen- 
tration at which the maximum in the conductivity is observed (9% and 13% for YSZ 
and ScSZ, respectively) and then two higher concentrations (18% and 24%). The rdfs 
for YSZ show that there is a preference for vacancies to be nearest neighbours to the 
small Zr'^^ cations, as found in the preceding paper for 11% YSZ. Interestingly, these 
binding tendencies become weaker as more dopant cations are added. This is probably 
caused by the fact that, as we add more Y2O3, the number of dopant cations increases 
(as well as the number of vacancies) and it becomes more difficult for the vacancies 
to avoid the dopant cations. In ScSZ, on the other hand, the vacancies have almost 
no preference for a certain cation species for all the dopant concentrations, as shown 
in figure |4j and this preference becomes even weaker as the dopant concentration is 
increased. Figure |4] seems therefore to confirm the hypothesis that cation- vacancy in- 
teractions alone cannot explain the anomalous behaviour of the conductivity in these 
materials. These effects are almost absent in ScSZ at the temperatures of interest and 
they become even less important at higher dopant concentrations. We now focus on the 
importance of vacancy-vacancy interactions. 

It might be anticipated that the vacancy-vacancy ordering would be affected by the 
vacancy-cation ordering. For this reason it is convenient to supplement the information 
on the vacancies in the real systems YSZ and ScSZ with that calculated on the ideal 
system i-YSZ;^, where all the cations are identical, so that the only vacancy-ordering 
effects are caused by the vacancy-vacancy interactions. The vacancy- vacancy rdfs for 
YSZ, ScSZ and i-YSZi have been studied as a function of the vacancy concentration 
and temperature. The peaks appear at the positions of the (simple cubic) oxide lattice 
sites, so that the first peak corresponds to a nearest-neighbour position along (100), the 
second to (110), the third to (111) etc. If the vacancies were randomly distributed over 
the lattice sites, with no correlations between them, the vacancy-vacancy rdf would 
have the same appearance as an anion-anion rdf (ga-a) for a perfect fluorite system 
(for instance PbF?) and we have included a plot of go_„ in both figures to illustrate 
the strength of the vacancy- vacancy correlations which are actually found: the relative 
intensities of the peaks in gv-v, relative to those in ga-a, give a guide to the probability 
of finding vacancies at these separations compared to the random distribution. 

The data for YSZ are illustrated in the top part of figure [5]for a simulation at 1250 
K and for three different dopant concentrations, x. At small separations it is clear that 
the (100) and (110) relative positions are very unfavourable relative to the random dis- 
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Figure 4; Cation-vacancy radial distribution functions for YSZ (top) and ScSZ (bottom) for 
different dopant concentrations at 1250 K. The red curves are the Zr-Vac rdfs while the black 
curves represent the dopant cation-vacancy rdfs. 



tribution, and that the occupancy of the (111) position is considerably enhanced. These 
propensities are consistent with the relative energy order which Bogicevic^ found in ab 
initio studies of the stability of different vacancy pairs at K. Indeed we have already 
seen in the preceding paper that the enthalpy difference between the <1 10> and <1 1 1> 
vacancy pairs coincides well with the energy difference found in the ab-initio studies. 
These propensities are also consistent with the interpretation^ of the diffuse scattering 
in the x = 0.24 sample discussed above. Also, these propensities do not seem to change 
as a function of the dopant concentrations studied here. It is interesting to note that 
occupancy of a pair of positions separated by the (111) vector results in a substantial 
lattice distortion, as seen in the position of the corresponding peak in gv-v compared 
to ga-a- It is also clear that the vacancy-vacancy correlations extend well beyond the 
first unit cell and, indeed, the rdf only begins to match that of the random distribution 
for separations larger than 9 A. The (200) position seems particularly unfavourable, but 
there is an enhanced occupation of the (210) and (211) positions which suggests that 
the pairs of vacancies are themselves beginning to order It is well known that at the 
composition Y/{Li:^0\2, corresponding to x = 0.4, the YSZ system forms a compound 
which can be described as based on the fluorite structure with ordered vacancy pairs in 
positions which are themselves ordered along the (112) direction and this appears to 
be the tendency which is being picked up in gv-v even at considerably smaller x values. 

The middle and bottom parts of figure |5] shows the vacancy-vacancy rdfs for ScSZ 
and i-YSZf . These show very similar trends to those observed in YSZ, with vacancies 
preferring pairing up in the (111) direction. A comparison between the top and middle 
part of figure [5] shows that these ordering tendencies are quite similar between the two 
materials and that therefore they are not influenced by the dopant cation species. The 
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implications of this finding are that the vacancy-vacancy ordering tendencies are not 
affected by the nature of the dopant species and that they are an intrinsic property of 
the fluorite lattice. 



5 Discussion and conclusions 

We have established that the strength of the vacancy-vacancy ordering effects is very 
similar in YSZ, ScSZ and i-YSZ^. Goff et al. showed that the vacancy clustering af- 
fects the conductivity of these materials by reducing the mobility of the vacancies'^ 
and we have confirmed this idea in the simulated systems. It remains to show at what 
dopant levels the vacancy-vacancy ordering effects begin to influence the conductivity. 

We can do this conveniently by calculating the jc-dependence of the conductivity 
of i-YSZ;c, in which there are no dopant-vacancy interactions, and comparing it with 
that of the real materials. This is shown in figure l6] The shape of the conductivity 
versus x curve for i-YSZv follows closely that of ScSZ, especially in the vicinity of 
the maximum, although the i-YSZ^ values are systematically higher At higher dopant 
concentrations, the two curves differ and this may reflect the fact that at these higher 
concentrations the equalisation of the charges on all the cations in i- YSZjc exerts a sig- 
nificant effect on the interactions and also the short-range potential in this system is, 
on average, slightly less repulsive than in the real systems. 

Figure |6] shows that, even if all the differences between the cation species are re- 
moved, a drop in the conductivity is still observed in i-YSZi and this happens at 13%, 
which is approximately the same vacancy concentration as in the "real" materials. This 
demonstrates that the observed anomalous behaviour is mainly caused by vacancy- 
vacancy interactions. Vacancy-vacancy interactions are caused by the fact that va- 
cancies are charged and distort the lattice. ^ As more vacancies are introduced in the 
material, they will try to minimise these interactions by ordering over the fluorite lat- 
tice. This limits the mobility of vacancies by correlating their motion, as explained in 
re f. 1211 Unlike cation-vacancy interactions, which can be reduced by choosing a suitable 
dopant cation which matches the radius of Zr'*+, vacancy- vacancy interactions can- 
not be minimised in such a way. This also implies that this is an intrinsic effect and 
therefore a "property" of the fluorite lattice. This explains the weak dependence of the 
maximum's position on the dopant cation's radius observed by Arachi et aZ. Dopants 
with larger radii (Y^+, Dy^+, Gd^+) will increase the cation-vacancy interaction and 
therefore slightly lower the maximum's position, but will not influence very much the 
vacancy-vacancy interactions, which is the main cause of the anomalous behaviour of 
the conductivity of these materials. For this reason ScSZ is probably the stabilised zir- 
conia with the highest achievable ionic conductivity, as the Sc^+ radius matches that of 
Zr^+ most closely. 
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Figure 5: Vacancy-vacancy radial distributiM functions for (Y20j,)x - (Zr02)i-x (top) and 
(Sc20:i)x - (Zr02)\-x (bottom) at 1250 K. 







xin(ZrO,),.^-(M,03)^ 



Figure 6: Conductivity versus vacancy concentration for i-YSZ^ , YSZ and ScSZ at 1670 K. 
The data for i-YSZx has been rescaled by a factor of 0.8 to facilitate the comparison. That i- 
YSZx has a slightly higher conductivity, at each x, than ScSZ is probably a consequence of the 
potential which has, on average, a less repulsive short-range interactions (see the discussion in 
reference^). 
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